home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / D-G / GemsII / inv_cmap / inv_cmap.c < prev   
Encoding:
C/C++ Source or Header  |  1992-06-16  |  19.5 KB  |  746 lines  |  [TEXT/MPS ]
  1. /*
  2.  * This software is copyrighted as noted below.  It may be freely copied,
  3.  * modified, and redistributed, provided that the copyright notice is
  4.  * preserved on all copies.
  5.  *
  6.  * There is no warranty or other guarantee of fitness for this software,
  7.  * it is provided solely "as is".  Bug reports or fixes may be sent
  8.  * to the author, who may or may not act on them as he desires.
  9.  *
  10.  * You may not include this software in a program or other software product
  11.  * without supplying the source, or without informing the end-user that the
  12.  * source is available for no extra charge.
  13.  *
  14.  * If you modify this software, you should include a notice giving the
  15.  * name of the person performing the modification, the date of modification,
  16.  * and the reason for such modification.
  17.  */
  18. /*
  19.  * inv_cmap.c - Compute an inverse colormap.
  20.  *
  21.  * Author:    Spencer W. Thomas
  22.  *         EECS Dept.
  23.  *         University of Michigan
  24.  * Date:    Thu Sep 20 1990
  25.  * Copyright (c) 1990, University of Michigan
  26.  *
  27.  * $Id: inv_cmap.c,v 3.0.1.1 90/11/29 15:09:43 spencer Exp $
  28.  */
  29.  
  30. #include <math.h>
  31. #include <stdio.h>
  32.  
  33. /* Print some performance stats. */
  34. /* #define INSTRUMENT_IT    */
  35. /* Track minimum and maximum in inv_cmap_2. */
  36. #define MINMAX_TRACK
  37.  
  38. static int bcenter, gcenter, rcenter;
  39. static long gdist, rdist, cdist;
  40. static long cbinc, cginc, crinc;
  41. static unsigned long *gdp, *rdp, *cdp;
  42. static unsigned char *grgbp, *rrgbp, *crgbp;
  43. static gstride, rstride;
  44. static long x, xsqr, colormax;
  45. static int cindex;
  46. #ifdef INSTRUMENT_IT
  47. static long outercount = 0, innercount = 0;
  48. #endif
  49.  
  50. /*****************************************************************
  51.  * TAG( inv_cmap_2 )
  52.  *
  53.  * Compute an inverse colormap efficiently.
  54.  * Inputs:
  55.  *     colors:        Number of colors in the forward colormap.
  56.  *     colormap:    The forward colormap.
  57.  *     bits:        Number of quantization bits.  The inverse
  58.  *             colormap will have (2^bits)^3 entries.
  59.  *     dist_buf:    An array of (2^bits)^3 long integers to be
  60.  *             used as scratch space.
  61.  * Outputs:
  62.  *     rgbmap:        The output inverse colormap.  The entry
  63.  *             rgbmap[(r<<(2*bits)) + (g<<bits) + b]
  64.  *             is the colormap entry that is closest to the
  65.  *             (quantized) color (r,g,b).
  66.  * Assumptions:
  67.  *     Quantization is performed by right shift (low order bits are
  68.  *     truncated).  Thus, the distance to a quantized color is
  69.  *     actually measured to the color at the center of the cell
  70.  *     (i.e., to r+.5, g+.5, b+.5, if (r,g,b) is a quantized color).
  71.  * Algorithm:
  72.  *     Uses a "distance buffer" algorithm:
  73.  *     The distance from each representative in the forward color map
  74.  *     to each point in the rgb space is computed.  If it is less
  75.  *     than the distance currently stored in dist_buf, then the
  76.  *     corresponding entry in rgbmap is replaced with the current
  77.  *     representative (and the dist_buf entry is replaced with the
  78.  *     new distance).
  79.  *
  80.  *     The distance computation uses an efficient incremental formulation.
  81.  *
  82.  *     Distances are computed "outward" from each color.  If the
  83.  *     colors are evenly distributed in color space, the expected
  84.  *     number of cells visited for color I is N^3/I.
  85.  *     Thus, the complexity of the algorithm is O(log(K) N^3),
  86.  *     where K = colors, and N = 2^bits.
  87.  */
  88.  
  89. /*
  90.  * Here's the idea:  scan from the "center" of each cell "out"
  91.  * until we hit the "edge" of the cell -- that is, the point
  92.  * at which some other color is closer -- and stop.  In 1-D,
  93.  * this is simple:
  94.  *     for i := here to max do
  95.  *         if closer then buffer[i] = this color
  96.  *         else break
  97.  *     repeat above loop with i := here-1 to min by -1
  98.  *
  99.  * In 2-D, it's trickier, because along a "scan-line", the
  100.  * region might start "after" the "center" point.  A picture
  101.  * might clarify:
  102.  *         |    ...
  103.  *               | ...    .
  104.  *              ...        .
  105.  *           ... |      .
  106.  *          .    +         .
  107.  *           .          .
  108.  *            .         .
  109.  *             .........
  110.  *
  111.  * The + marks the "center" of the above region.  On the top 2
  112.  * lines, the region "begins" to the right of the "center".
  113.  *
  114.  * Thus, we need a loop like this:
  115.  *     detect := false
  116.  *     for i := here to max do
  117.  *         if closer then
  118.  *             buffer[..., i] := this color
  119.  *             if !detect then
  120.  *                 here = i
  121.  *                 detect = true
  122.  *         else
  123.  *             if detect then
  124.  *                 break
  125.  *                 
  126.  * Repeat the above loop with i := here-1 to min by -1.  Note that
  127.  * the "detect" value should not be reinitialized.  If it was
  128.  * "true", and center is not inside the cell, then none of the
  129.  * cell lies to the left and this loop should exit
  130.  * immediately.
  131.  *
  132.  * The outer loops are similar, except that the "closer" test
  133.  * is replaced by a call to the "next in" loop; its "detect"
  134.  * value serves as the test.  (No assignment to the buffer is
  135.  * done, either.)
  136.  *
  137.  * Each time an outer loop starts, the "here", "min", and
  138.  * "max" values of the next inner loop should be
  139.  * re-initialized to the center of the cell, 0, and cube size,
  140.  * respectively.  Otherwise, these values will carry over from
  141.  * one "call" to the inner loop to the next.  This tracks the
  142.  * edges of the cell and minimizes the number of
  143.  * "unproductive" comparisons that must be made.
  144.  *
  145.  * Finally, the inner-most loop can have the "if !detect"
  146.  * optimized out of it by splitting it into two loops: one
  147.  * that finds the first color value on the scan line that is
  148.  * in this cell, and a second that fills the cell until
  149.  * another one is closer:
  150.  *      if !detect then        {needed for "down" loop}
  151.  *         for i := here to max do
  152.  *         if closer then
  153.  *             buffer[..., i] := this color
  154.  *             detect := true
  155.  *             break
  156.  *    for i := i+1 to max do
  157.  *        if closer then
  158.  *             buffer[..., i] := this color
  159.  *         else
  160.  *             break
  161.  *
  162.  * In this implementation, each level will require the
  163.  * following variables.  Variables labelled (l) are local to each
  164.  * procedure.  The ? should be replaced with r, g, or b:
  165.  *      cdist:            The distance at the starting point.
  166.  *     ?center:    The value of this component of the color
  167.  *      c?inc:            The initial increment at the ?center position.
  168.  *     ?stride:    The amount to add to the buffer
  169.  *             pointers (dp and rgbp) to get to the
  170.  *             "next row".
  171.  *     min(l):        The "low edge" of the cell, init to 0
  172.  *     max(l):        The "high edge" of the cell, init to
  173.  *             colormax-1
  174.  *     detect(l):        True if this row has changed some
  175.  *                     buffer entries.
  176.  *      i(l):             The index for this row.
  177.  *      ?xx:            The accumulated increment value.
  178.  *      
  179.  *      here(l):        The starting index for this color.  The
  180.  *                      following variables are associated with here,
  181.  *                      in the sense that they must be updated if here
  182.  *                      is changed.
  183.  *      ?dist:            The current distance for this level.  The
  184.  *                      value of dist from the previous level (g or r,
  185.  *                      for level b or g) initializes dist on this
  186.  *                      level.  Thus gdist is associated with here(b)).
  187.  *      ?inc:            The initial increment for the row.
  188.  
  189.  *      ?dp:            Pointer into the distance buffer.  The value
  190.  *                      from the previous level initializes this level.
  191.  *      ?rgbp:            Pointer into the rgb buffer.  The value
  192.  *                      from the previous level initializes this level.
  193.  * 
  194.  * The blue and green levels modify 'here-associated' variables (dp,
  195.  * rgbp, dist) on the green and red levels, respectively, when here is
  196.  * changed.
  197.  */
  198.  
  199. void
  200. inv_cmap_2( colors, colormap, bits, dist_buf, rgbmap )
  201. int colors, bits;
  202. unsigned char *colormap[3], *rgbmap;
  203. unsigned long *dist_buf;
  204. {
  205.     int nbits = 8 - bits;
  206.  
  207.     colormax = 1 << bits;
  208.     x = 1 << nbits;
  209.     xsqr = 1 << (2 * nbits);
  210.  
  211.     /* Compute "strides" for accessing the arrays. */
  212.     gstride = colormax;
  213.     rstride = colormax * colormax;
  214.  
  215. #ifdef INSTRUMENT_IT
  216.     outercount = 0;
  217.     innercount = 0;
  218. #endif
  219.     maxfill( dist_buf, colormax );
  220.  
  221.     for ( cindex = 0; cindex < colors; cindex++ )
  222.     {
  223.     /*
  224.      * Distance formula is
  225.      * (red - map[0])^2 + (green - map[1])^2 + (blue - map[2])^2
  226.      *
  227.      * Because of quantization, we will measure from the center of
  228.      * each quantized "cube", so blue distance is
  229.      *     (blue + x/2 - map[2])^2,
  230.      * where x = 2^(8 - bits).
  231.      * The step size is x, so the blue increment is
  232.      *     2*x*blue - 2*x*map[2] + 2*x^2
  233.      *
  234.      * Now, b in the code below is actually blue/x, so our
  235.      * increment will be 2*(b*x^2 + x^2 - x*map[2]).  For
  236.      * efficiency, we will maintain this quantity in a separate variable
  237.      * that will be updated incrementally by adding 2*x^2 each time.
  238.      */
  239.     /* The initial position is the cell containing the colormap
  240.      * entry.  We get this by quantizing the colormap values.
  241.      */
  242.     rcenter = colormap[0][cindex] >> nbits;
  243.     gcenter = colormap[1][cindex] >> nbits;
  244.     bcenter = colormap[2][cindex] >> nbits;
  245.  
  246.     rdist = colormap[0][cindex] - (rcenter * x + x/2);
  247.     gdist = colormap[1][cindex] - (gcenter * x + x/2);
  248.     cdist = colormap[2][cindex] - (bcenter * x + x/2);
  249.     cdist = rdist*rdist + gdist*gdist + cdist*cdist;
  250.  
  251.     crinc = 2 * ((rcenter + 1) * xsqr - (colormap[0][cindex] * x));
  252.     cginc = 2 * ((gcenter + 1) * xsqr - (colormap[1][cindex] * x));
  253.     cbinc = 2 * ((bcenter + 1) * xsqr - (colormap[2][cindex] * x));
  254.  
  255.     /* Array starting points. */
  256.     cdp = dist_buf + rcenter * rstride + gcenter * gstride + bcenter;
  257.     crgbp = rgbmap + rcenter * rstride + gcenter * gstride + bcenter;
  258.  
  259.     (void)redloop();
  260.     }
  261. #ifdef INSTRUMENT_IT
  262.     fprintf( stderr, "K = %d, N = %d, outer count = %ld, inner count = %ld\n",
  263.          colors, colormax, outercount, innercount );
  264. #endif
  265. }
  266.  
  267. /* redloop -- loop up and down from red center. */
  268. int
  269. redloop()
  270. {
  271.     int detect;
  272.     int r, i = cindex;
  273.     int first;
  274.     long txsqr = xsqr + xsqr;
  275.     static int here, min, max;
  276.     static long rxx;
  277.  
  278.     detect = 0;
  279.  
  280.     /* Basic loop up. */
  281.     for ( r = rcenter, rdist = cdist, rxx = crinc,
  282.       rdp = cdp, rrgbp = crgbp, first = 1;
  283.       r < colormax;
  284.       r++, rdp += rstride, rrgbp += rstride,
  285.       rdist += rxx, rxx += txsqr, first = 0 )
  286.     {
  287.     if ( greenloop( first ) )
  288.         detect = 1;
  289.     else if ( detect )
  290.         break;
  291.     }
  292.     
  293.     /* Basic loop down. */
  294.     for ( r = rcenter - 1, rxx = crinc - txsqr, rdist = cdist - rxx,
  295.       rdp = cdp - rstride, rrgbp = crgbp - rstride, first = 1;
  296.       r >= 0;
  297.       r--, rdp -= rstride, rrgbp -= rstride,
  298.       rxx -= txsqr, rdist -= rxx, first = 0 )
  299.     {
  300.     if ( greenloop( first ) )
  301.         detect = 1;
  302.     else if ( detect )
  303.         break;
  304.     }
  305.     
  306.     return detect;
  307. }
  308.  
  309. /* greenloop -- loop up and down from green center. */
  310. int
  311. greenloop( restart )
  312. {
  313.     int detect;
  314.     int g, i = cindex;
  315.     int first;
  316.     long txsqr = xsqr + xsqr;
  317.     static int here, min, max;
  318. #ifdef MINMAX_TRACK
  319.     static int prevmax, prevmin;
  320.     int thismax, thismin;
  321. #endif
  322.     static long ginc, gxx, gcdist;    /* "gc" variables maintain correct */
  323.     static unsigned long *gcdp;        /*  values for bcenter position, */
  324.     static unsigned char *gcrgbp;    /*  despite modifications by blueloop */
  325.                     /*  to gdist, gdp, grgbp. */
  326.  
  327.     if ( restart )
  328.     {
  329.     here = gcenter;
  330.     min = 0;
  331.     max = colormax - 1;
  332.     ginc = cginc;
  333. #ifdef MINMAX_TRACK
  334.     prevmax = 0;
  335.     prevmin = colormax;
  336. #endif
  337.     }
  338.  
  339. #ifdef MINMAX_TRACK
  340.     thismin = min;
  341.     thismax = max;
  342. #endif
  343.     detect = 0;
  344.  
  345.     /* Basic loop up. */
  346.     for ( g = here, gcdist = gdist = rdist, gxx = ginc,
  347.       gcdp = gdp = rdp, gcrgbp = grgbp = rrgbp, first = 1;
  348.       g <= max;
  349.       g++, gdp += gstride, gcdp += gstride, grgbp += gstride, gcrgbp += gstride,
  350.       gdist += gxx, gcdist += gxx, gxx += txsqr, first = 0 )
  351.     {
  352.     if ( blueloop( first ) )
  353.     {
  354.         if ( !detect )
  355.         {
  356.         /* Remember here and associated data! */
  357.         if ( g > here )
  358.         {
  359.             here = g;
  360.             rdp = gcdp;
  361.             rrgbp = gcrgbp;
  362.             rdist = gcdist;
  363.             ginc = gxx;
  364. #ifdef MINMAX_TRACK
  365.             thismin = here;
  366. #endif
  367.         }
  368.         detect = 1;
  369.         }
  370.     }
  371.     else if ( detect )
  372.     {
  373. #ifdef MINMAX_TRACK
  374.         thismax = g - 1;
  375. #endif
  376.         break;
  377.     }
  378.     }
  379.     
  380.     /* Basic loop down. */
  381.     for ( g = here - 1, gxx = ginc - txsqr, gcdist = gdist = rdist - gxx,
  382.       gcdp = gdp = rdp - gstride, gcrgbp = grgbp = rrgbp - gstride,
  383.       first = 1;
  384.       g >= min;
  385.       g--, gdp -= gstride, gcdp -= gstride, grgbp -= gstride, gcrgbp -= gstride,
  386.       gxx -= txsqr, gdist -= gxx, gcdist -= gxx, first = 0 )
  387.     {
  388.     if ( blueloop( first ) )
  389.     {
  390.         if ( !detect )
  391.         {
  392.         /* Remember here! */
  393.         here = g;
  394.         rdp = gcdp;
  395.         rrgbp = gcrgbp;
  396.         rdist = gcdist;
  397.         ginc = gxx;
  398. #ifdef MINMAX_TRACK
  399.         thismax = here;
  400. #endif
  401.         detect = 1;
  402.         }
  403.     }
  404.     else if ( detect )
  405.     {
  406. #ifdef MINMAX_TRACK
  407.         thismin = g + 1;
  408. #endif
  409.         break;
  410.     }
  411.     }
  412.     
  413. #ifdef MINMAX_TRACK
  414.     /* If we saw something, update the edge trackers.  For now, only
  415.      * tracks edges that are "shrinking" (min increasing, max
  416.      * decreasing.
  417.      */
  418.     if ( detect )
  419.     {
  420.     if ( thismax < prevmax )
  421.         max = thismax;
  422.  
  423.     prevmax = thismax;
  424.  
  425.     if ( thismin > prevmin )
  426.         min = thismin;
  427.  
  428.     prevmin = thismin;
  429.     }
  430. #endif
  431.  
  432.     return detect;
  433. }
  434.  
  435. /* blueloop -- loop up and down from blue center. */
  436. int
  437. blueloop( restart )
  438. {
  439.     int detect;
  440.     register unsigned long *dp;
  441.     register unsigned char *rgbp;
  442.     register long bdist, bxx;
  443.     register int b, i = cindex;
  444.     register long txsqr = xsqr + xsqr;
  445.     register int lim;
  446.     static int here, min, max;
  447. #ifdef MINMAX_TRACK
  448.     static int prevmin, prevmax;
  449.     int thismin, thismax;
  450. #endif /* MINMAX_TRACK */
  451.     static long binc;
  452.  
  453.     if ( restart )
  454.     {
  455.     here = bcenter;
  456.     min = 0;
  457.     max = colormax - 1;
  458.     binc = cbinc;
  459. #ifdef MINMAX_TRACK
  460.     prevmin = colormax;
  461.     prevmax = 0;
  462. #endif /* MINMAX_TRACK */
  463.     }
  464.  
  465.     detect = 0;
  466. #ifdef MINMAX_TRACK
  467.     thismin = min;
  468.     thismax = max;
  469. #endif
  470.  
  471.     /* Basic loop up. */
  472.     /* First loop just finds first applicable cell. */
  473.     for ( b = here, bdist = gdist, bxx = binc, dp = gdp, rgbp = grgbp, lim = max;
  474.       b <= lim;
  475.       b++, dp++, rgbp++,
  476.       bdist += bxx, bxx += txsqr )
  477.     {
  478. #ifdef INSTRUMENT_IT
  479.     outercount++;
  480. #endif
  481.     if ( *dp > bdist )
  482.     {
  483.         /* Remember new 'here' and associated data! */
  484.         if ( b > here )
  485.         {
  486.         here = b;
  487.         gdp = dp;
  488.         grgbp = rgbp;
  489.         gdist = bdist;
  490.         binc = bxx;
  491. #ifdef MINMAX_TRACK
  492.         thismin = here;
  493. #endif
  494.         }
  495.         detect = 1;
  496. #ifdef INSTRUMENT_IT
  497.         outercount--;
  498. #endif
  499.         break;
  500.     }
  501.     }
  502.     /* Second loop fills in a run of closer cells. */
  503.     for ( ;
  504.       b <= lim;
  505.       b++, dp++, rgbp++,
  506.       bdist += bxx, bxx += txsqr )
  507.     {
  508. #ifdef INSTRUMENT_IT
  509.     outercount++;
  510. #endif
  511.     if ( *dp > bdist )
  512.     {
  513. #ifdef INSTRUMENT_IT
  514.         innercount++;
  515. #endif
  516.         *dp = bdist;
  517.         *rgbp = i;
  518.     }
  519.     else
  520.     {
  521. #ifdef MINMAX_TRACK
  522.         thismax = b - 1;
  523. #endif
  524.         break;
  525.     }
  526.     }
  527.     
  528.     /* Basic loop down. */
  529.     /* Do initializations here, since the 'find' loop might not get
  530.      * executed. 
  531.      */
  532.     lim = min;
  533.     b = here - 1;
  534.     bxx = binc - txsqr;
  535.     bdist = gdist - bxx;
  536.     dp = gdp - 1;
  537.     rgbp = grgbp - 1;
  538.     /* The 'find' loop is executed only if we didn't already find
  539.      * something.
  540.      */
  541.     if ( !detect )
  542.     for ( ;
  543.           b >= lim;
  544.           b--, dp--, rgbp--,
  545.           bxx -= txsqr, bdist -= bxx )
  546.     {
  547. #ifdef INSTRUMENT_IT
  548.         outercount++;
  549. #endif
  550.         if ( *dp > bdist )
  551.         {
  552.         /* Remember here! */
  553.         /* No test for b against here necessary because b <
  554.          * here by definition.
  555.          */
  556.         here = b;
  557.         gdp = dp;
  558.         grgbp = rgbp;
  559.         gdist = bdist;
  560.         binc = bxx;
  561. #ifdef MINMAX_TRACK
  562.         thismax = here;
  563. #endif
  564.         detect = 1;
  565. #ifdef INSTRUMENT_IT
  566.         outercount--;
  567. #endif
  568.         break;
  569.         }
  570.     }
  571.     /* The 'update' loop. */
  572.     for ( ;
  573.       b >= lim;
  574.       b--, dp--, rgbp--,
  575.       bxx -= txsqr, bdist -= bxx )
  576.     {
  577. #ifdef INSTRUMENT_IT
  578.     outercount++;
  579. #endif
  580.     if ( *dp > bdist )
  581.     {
  582. #ifdef INSTRUMENT_IT
  583.         innercount++;
  584. #endif
  585.         *dp = bdist;
  586.         *rgbp = i;
  587.     }
  588.     else
  589.     {
  590. #ifdef MINMAX_TRACK
  591.         thismin = b + 1;
  592. #endif
  593.         break;
  594.     }
  595.     }
  596.  
  597.  
  598.     /* If we saw something, update the edge trackers. */
  599. #ifdef MINMAX_TRACK
  600.     if ( detect )
  601.     {
  602.     /* Only tracks edges that are "shrinking" (min increasing, max
  603.      * decreasing.
  604.      */
  605.     if ( thismax < prevmax )
  606.         max = thismax;
  607.  
  608.     if ( thismin > prevmin )
  609.         min = thismin;
  610.     
  611.     /* Remember the min and max values. */
  612.     prevmax = thismax;
  613.     prevmin = thismin;
  614.     }
  615. #endif /* MINMAX_TRACK */
  616.  
  617.     return detect;
  618. }
  619.  
  620. maxfill( buffer, side )
  621. unsigned long *buffer;
  622. long side;
  623. {
  624.     register unsigned long maxv = ~0L;
  625.     register long i;
  626.     register unsigned long *bp;
  627.  
  628.     for ( i = colormax * colormax * colormax, bp = buffer;
  629.       i > 0;
  630.       i--, bp++ )
  631.     *bp = maxv;
  632. }
  633.  
  634. /*****************************************************************
  635.  * TAG( inv_cmap_1 )
  636.  *
  637.  * Compute an inverse colormap efficiently.
  638.  * Inputs:
  639.  
  640.  *     colors:        Number of colors in the forward colormap.
  641.  *     colormap:    The forward colormap.
  642.  *     bits:        Number of quantization bits.  The inverse
  643.  *             colormap will have (2^bits)^3 entries.
  644.  *     dist_buf:    An array of (2^bits)^3 long integers to be
  645.  *             used as scratch space.
  646.  * Outputs:
  647.  *     rgbmap:        The output inverse colormap.  The entry
  648.  *             rgbmap[(r<<(2*bits)) + (g<<bits) + b]
  649.  *             is the colormap entry that is closest to the
  650.  *             (quantized) color (r,g,b).
  651.  * Assumptions:
  652.  *     Quantization is performed by right shift (low order bits are
  653.  *     truncated).  Thus, the distance to a quantized color is
  654.  *     actually measured to the color at the center of the cell
  655.  *     (i.e., to r+.5, g+.5, b+.5, if (r,g,b) is a quantized color).
  656.  * Algorithm:
  657.  *     Uses a "distance buffer" algorithm:
  658.  *     The distance from each representative in the forward color map
  659.  *     to each point in the rgb space is computed.  If it is less
  660.  *     than the distance currently stored in dist_buf, then the
  661.  *     corresponding entry in rgbmap is replaced with the current
  662.  *     representative (and the dist_buf entry is replaced with the
  663.  *     new distance).
  664.  *
  665.  *     The distance computation uses an efficient incremental formulation.
  666.  *
  667.  *     Right now, distances are computed for all entries in the rgb
  668.  *     space.  Thus, the complexity of the algorithm is O(K N^3),
  669.  *     where K = colors, and N = 2^bits.
  670.  */
  671. void
  672. inv_cmap_1( colors, colormap, bits, dist_buf, rgbmap )
  673. int colors, bits;
  674. unsigned char *colormap[3], *rgbmap;
  675. unsigned long *dist_buf;
  676. {
  677.     register unsigned long *dp;
  678.     register unsigned char *rgbp;
  679.     register long bdist, bxx;
  680.     register int b, i;
  681.     int nbits = 8 - bits;
  682.     register int colormax = 1 << bits;
  683.     register long xsqr = 1 << (2 * nbits);
  684.     int x = 1 << nbits;
  685.     int rinc, ginc, binc, r, g;
  686.     long rdist, gdist, rxx, gxx;
  687.  
  688.     for ( i = 0; i < colors; i++ )
  689.     {
  690.     /*
  691.      * Distance formula is
  692.      * (red - map[0])^2 + (green - map[1])^2 + (blue - map[2])^2
  693.      *
  694.      * Because of quantization, we will measure from the center of
  695.      * each quantized "cube", so blue distance is
  696.      *     (blue + x/2 - map[2])^2,
  697.      * where x = 2^(8 - bits).
  698.      * The step size is x, so the blue increment is
  699.      *     2*x*blue - 2*x*map[2] + 2*x^2
  700.      *
  701.      * Now, b in the code below is actually blue/x, so our
  702.      * increment will be 2*x*x*b + (2*x^2 - 2*x*map[2]).  For
  703.      * efficiency, we will maintain this quantity in a separate variable
  704.      * that will be updated incrementally by adding 2*x^2 each time.
  705.      */
  706.     rdist = colormap[0][i] - x/2;
  707.     gdist = colormap[1][i] - x/2;
  708.     bdist = colormap[2][i] - x/2;
  709.     rdist = rdist*rdist + gdist*gdist + bdist*bdist;
  710.  
  711.     rinc = 2 * (xsqr - (colormap[0][i] << nbits));
  712.     ginc = 2 * (xsqr - (colormap[1][i] << nbits));
  713.     binc = 2 * (xsqr - (colormap[2][i] << nbits));
  714.     dp = dist_buf;
  715.     rgbp = rgbmap;
  716.     for ( r = 0, rxx = rinc;
  717.           r < colormax;
  718.           rdist += rxx, r++, rxx += xsqr + xsqr )
  719.         for ( g = 0, gdist = rdist, gxx = ginc;
  720.           g < colormax;
  721.           gdist += gxx, g++, gxx += xsqr + xsqr )
  722.         for ( b = 0, bdist = gdist, bxx = binc;
  723.               b < colormax;
  724.               bdist += bxx, b++, dp++, rgbp++,
  725.               bxx += xsqr + xsqr )
  726.         {
  727. #ifdef INSTRUMENT_IT
  728.             outercount++;
  729. #endif
  730.             if ( i == 0 || *dp > bdist )
  731.             {
  732. #ifdef INSTRUMENT_IT
  733.             innercount++;
  734. #endif
  735.             *dp = bdist;
  736.             *rgbp = i;
  737.             }
  738.         }
  739.     }
  740. #ifdef INSTRUMENT_IT
  741.     fprintf( stderr, "K = %d, N = %d, outer count = %ld, inner count = %ld\n",
  742.          colors, colormax, outercount, innercount );
  743. #endif
  744. }
  745.  
  746.